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^A  new  approach  for  the  generation  of  flow-adaptive 
grids  for  numerical  solution  of  fluid  dynamics  problems  is 
presented.  However,  the  method  is  applicable  to  the  numer¬ 
ical  evaluation  of  any  partial  differential  equation. 

The  dynamic  coupling  of  the  grid  with  the  flow  solu 
tion  is  accomplished  through  a  grid-optimization  technique. 
The  optimization  is  based  on  the  minimization  of  the  finite 
difference  truncation  error  in  the  transformed  plane.  The 
method  is  tested  on  the  one-dimensional  Burgers'  equation 

which  is  representative  of  typical  fluid  dynamics  problems. 

0«<*  ■*- 

Burgers'  equation  is  solved  with  an  optimized  SOR^method 
using  upwind  differences  for  the  convective  term. 

Results  are  presented  for  various  Reynolds  numbers 
and  are  compared  to  results  from  a  similar  adaptive  grid 
method  and  to  results  for  a  static  grid.  They  show  the 
ability  of  the  method  to  concentrate  grid  points  high  in 
gradient  regions  where  large  truncation  errors  occur. 
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ADAPTIVE  GRID  GENERATION  FOR  NUMERICAL  SOLUTION 
OF  PARTIAL  DIFFERENTIAL  EQUATIONS 

I.  INTRODUCTION 

Fluid  mechanics  and  heat  transfer  problems  are  charac¬ 
terized  by  complex  nonlinear  partial  differential  equations, 
for  which  analytic  solutions  can  be  obtained  for  only  a  few 
limited  cases.  Due  to  the  rapid  increase  in  speed  and  mem¬ 
ory,  and  declining  cost  of  the  digital  computer,  an  ever  in¬ 
creasing  emphasis  is  being  placed  on  numerical  solution  of 
the  governing  differential  equations  by  finite  difference 
methods.  In  the  past  two  decades,  considerable  progress  has 
been  made  on  the  development  of  more  efficient  and  accurate 
numerical  algorithms. 

A  vital  component  of  any  numerical  algorithm  is  the  grid 
upon  which  the  solution  is  obtained.  The  solution  can  be 
greatly  simplified  if  a  well-suited  grid  is  chosen,  just  as 
the  choice  of  cylindrical  coordinates  rather  than  rectangu¬ 
lar  coordinates  simplifies  the  solution  of  the  potential  equa 
tions  about  a  cylinder.  In  contrast,  a  poorly  constructed 
grid  can  not  only  lead  to  large  errors  in  the  solution,  but 
also  can  cause  instability  resulting  in  a  totally  incorrect 
solution  or  no  solution  at  all.  The  area  of  numerical  grid 
generation  is  relatively  young  receiving  much  attention  in 
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the  past  ten  years.  As  J.  F.  Thompson,  a  leader  in  the  field, 
writes  (1:1) 

This  area  involves  the  engineer's  feel  for  the 
physical  behavior,  the  mathematician's  under¬ 
standing  of  the  functional  behavior,  and  a  lot 
of  imagination,  with  perhaps  a  little  help  from 
Urania. 

Two  very  important  factors  play  a  role  in  the  choice  of 
the  grid.  First,  the  coordinates  should  be  .surface  oriented 
so  that  the  surface  boundary  conditions  may  be  implemented 
with  no  need  for  interpolation  between  grid  points.  This 
not  only  simplifies  the  boundary  condition  implementation 
but  also  eliminates  large  errors  that  may  be  produced  in  the 
interpolation  process.  In  addition,  surface  oriented  coor¬ 
dinates  also  permit  coordinate-related  approximations  in  the 
flow  equations  (2:35).  Second,  the  grid  must  accurately  re¬ 
solve  high  gradient  regions  which  are  common  in  fluid  dynam¬ 
ics  problems,  because  it  is  in  these  regions  where  large 
numerical  errors  occur. 

One  solution  of  the  first  problem  is  the  use  of  partial 
differential  equations  to  generate  the  grid.  This  method 
was  popularized  by  Thompson,  Thames,  and  Mas  tin  in  1974  (3) 
and  is  widely  used  today.  The  grid  generation  procedure  in 
this  thesis,  presented  in  more  detail  later,  is  based  on  a 
modified  version  of  the  original  equations  presented  in  that 
paper . 

The  solution  of  the  second  problem  is  more  difficult. 
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The  method  of  Thompson  et  al.  provides  control  for  the  stretch- 
’  ing  of  interior  regions  to  resolve  the  high  gradients;  how¬ 

ever,  there  is,  in  general,  no  method  for  determining  the  val¬ 
ues  of  those  control  functions  due  to  the  lack  of  a  priori  in¬ 
formation  concerning  those  gradients  (4:28).  In  Ref.  4,  a 
method  was  developed  to  compute  the  forcing  function  to  mini¬ 
mize  truncation  error  using  boundary  layer  theory;  the  result¬ 
ing  mesh  system  is  referred  to  as  a  boundary  layer  dependent 
coordinate  system.  Therefore,  it  is  only  valid  in  regions  in 
which  boundary  layer  theory  is  applicable. 

In  recent  years,  flow  adaptive  grids  in  which  the  grid 
point  distribution  is  dynamically  coupled  to  the  developing 
solution  have  emerged  as  a  means  to  tackle  the.  second  problem. 
Several  different  methods  have  been  developed  to  achieve  the 
same  ef fect--concentrate  coordinates  in  high  gradient  regions 
thereby  decreasing  truncation  error  and  minimizing  the  number 
of  grid  points  necessary  to  produce  a  satisfactory  solution. 
Pierson  and  Kutler  (5)  describe  a  method  in  which  the  grid  is 
defined  by  a  minimization  of  the  local  truncation  error  in 
the  Least  squares  sense.  The  grid  is  then  algebraically  gen¬ 
erated,  using  Chebyshev  polynomials.  Saltzman  and  Brackbill 
(6)  define  their  grid  based  on  a  variational  analysis  result¬ 
ing  in  a  system  of  partial  differential  equations  whose  solu¬ 
tion  produces  the  grid.  In  the  work  of  Dwyer  et  al.  (7)  the 
grid  points  were  moved  in  time,  based  on  the  gradients  in  the 
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flow  variabLe.  In  the  paper  by  Ghia  et  al.  (2),  the  grid  a- 
daption  criterion  is  based  on  the  minimization  of  the  coef¬ 
ficient  of  the  convective  term  in  the  transformed  flow  equa¬ 
tions.  Freeman  (8)  describes  a  method  used  in  conjunction  with 
the  Thompson  elliptic  grid  generation  equations  in  which  the 
grid  control  functions  are  determined,  based  on  solution  grad¬ 
ients.  Anderson  and  Rai  (9)  describe  another  method  in  which 
the  grid  points  move  directly  under  an  attractive/repulsive  in¬ 
fluence  of  one  another,  based  on  the  magnitude  of  the  local  er¬ 
ror  compared  to  the  global  error.  This  influence  can  be  com¬ 
pared  to  the  force  of  a  test  charge  in  an  electrostatic  field 
or  to  a  gravitational  field  which  can  repel  as  well  as  attract. 

Anderson  and  Rai  list  the  following  considerations  for  the 
development  of  an  adaptive  grid  (9:320): 

1.  The  grid  must  evolve  as  part  of  the  solution. 

2.  Grid  points  must  move  due  to  both  boundary  motion  and 
changes  in  the  interior  solution. 

3.  The  grid  speed  equations  should  be  as  simple  as 
possible . 

A.  The  grid  speed  equations  must  account  for  the  ellip¬ 
tic  nature  of  the  problem. 

5.  The  resulting  grid  must  reduce  error,  provide  better 
resolution,  or  otherwise  improve  the  solution. 

6.  The  adaptive  grid  scheme  must  be  easily  extended  to 
any  number  of  dimensions. 
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In  addition,  it  is  felt  that  the  algorithm  should  be  robust 
enough  to  handle  a  variety  of  problems  with  arbitrary  input 
data . 

The  objective  of  this  thesis  is  to  develop  an  adaptive 
grid  which  is  based  on  the  systematic  determination  of  the 
grid  generation  control  functions  used  by  the  elliptic  grid 
generation  equations  of  Thompson  et  al.  Another  goal  of  this 
study  is  that,  in  the  process  of  developing  the  adaptive  grid, 
an  "optimum  grid"  is  produced  which  reduces  truncation  error, 
thus  satisfying  all  conditions  of  item  5  above. 
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I I .  Mathematical  Formulation 

The  adaptive  grid  method  presented  in  this  thesis  is 
based  on  a  one-dimensional  analysis  for  ease  of  formulation 
and  cost  considerations.  However,  par^s  of  the  following 
analysis  include  the  second  dimension  which  is  necessary  to 
present  a  clear  background  for  the  method  developed.  In  the 
following,  subscripts  denote  partial  differentiation. 

Background 

For  complex  geometries,  it  is  convenient  to  transform 
the  governing  differential  equations  from  the  physical  plane 
to  a  computational  plane  where  the  solution  is  obtained  and 
then  transformed  back  to  the  physical  plane.  For  one-dimen¬ 
sional,  time  dependent  problems,  this  transformation  is 

t  =  t  C  -  C(x,  t)  (1) 

where  x,  t  are  the  physical  variables  and  £,t  are  the  com¬ 
putational  variables.  The  derivatives  in  the  physical  plane 
transform  as 


=  f5/J 

(2a) 

=  r  f  _  _Aj  f  1/ 

x?  r5J/J 

(2b) 

=  fT  -  fc  XT  /J 

(2c) 
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where  if  is  the  dependent  flow  variable  and  J  is  the  Jacob¬ 
ian  of  the  transformation  given  by 

J  *  X£  =  l/€x  (2d) 

Terms  involving  derivatives  of  physical  space  coordinates 
with  respect  to  computational  space  coordinates  or  vice  versa 
are  referred  to  as  the  metrics  of  the  transformation  and  x 

T 

is  referred  to  as  the  grid  speed.  Due  to  the  metrics,  this 
transformation  renders  the  governing  equations  quite  complex; 
however,  there  are  three  attractive  reasons  which  out  weigh 
the  added  complexity.  First,  in  the  computational  plane,  the 
solution  can  be  performed  on  a  fixed  rectangular  grid  with 
uniform  spacing.  Second,  the  transformation  makes  it  possible 
to  concentrate  the  distribution  of  curvilinear  grid  lines  in 
the  physical  plane  in  regions  of  high  radients.  Third,  and 
most  importantly,  the  grid  lines  can  be  made  to  correspond  to 
the  boundaries  in  the  physical  plane,  no  matter  what  the  shape 
Figure  1  shows  the  idea  of  the  general  coordinate  transforma¬ 
tion  . 

Up  to  now,  the  actual  transformation  has  yet  to  be  speci¬ 
fied.  One  of  the  most  popular  methods  for  defining  this  trans 
formation  is  the  method  of  Thompson  et  al.  (3).  According  to 
this  method,  the  grid  is  determined  as  the  solution  of  a  set 
of  elliptic  differential  equations 


xx 


yy 


P(5,n) 


(3a) 
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'’xx  +  V  =  Q(i”')  (3b) 

where  x,  y  are  the  physical  corrdinates,  n  are  the  compu¬ 
tational  coordinates,  and  P,  Q  are  functions  which  control 
the  grid  spacing  in  the  interior  of  the  region,  hereafter 
called  grid  control  functions. 

The  solution  of  Eqs  (3)  may  be  no  easier  to  obtain  than 
that  of  the  flow  equations.  However,  if  the  roles  of  the  de¬ 
pendent  and  independent  variables  are  interchanged  so  that 
the  solution  is  performed  in  the  computational  plane,  the 
boundary  conditions  may  be  specified  along  constant  values  of 
the  computational  coordinates.  This  results  in 


“x«  •  2Bxen  +  yxnn 


-Jx(PxE  +  Qy^) 


^EE  -  26yEr,  +  yyn-,  ’  'J  (PxE  +  Qyn) 


where 


2  2 
a  =  x  +  y 

n  n 


8  =  xs  Xn  +  h  yn 


V  2  X  2 

Y  =  x  +  y 

S  £ 


y?  yn  +  yn  ^ 


(4a) 

(4b) 

(4c) 

(  4d) 

(4e) 

(  4f  ) 
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The  Jacobian  J  must  be  non-zero  for  a  unique  transformation. 
For  one  dimension,  Eqs  (4)  reduce  to 

xc;  "  Px53  '  0  (5) 

Equations  (4)  and  (5)  are  highly  nonlinear;  their  solu¬ 
tion  may  lead  to  numerical  difficulties  in  terms  of  oscilla¬ 
tions  and  instabilities  for  even  moderate  values  of  the  con¬ 
trol  functions  P,  Q  (4:21).  An  alernate  set  of  grid  genera¬ 
tion  equations  has  been  proposed  (11:57)  which  eliminates 
this  nonlinearity.  The  new  equations  are 


fr  l 

X 

X 

+  £ 

yy 

<5x2  + 

Cy2>  p 

(€  » n) 

(6a) 

n 

'  XX 

+  nyy 

<nx2  + 

Hy2)  Q 

U  ,  n) 

(6b) 

which, 

after  inversion, 

give 

ax 

20x  + 

Sn 

yx  * 

nn 

-  (aPx^ 

+  yQx 

T 

) 

i 

(7a) 

ay55  - 

2fJy5n  + 

yy 

nn 

-  (aPy^ 

+  yQy  ^ 

) 

i 

(7b) 

where  a,  8,  y  are  the  same  as  in  Eqs  (4). 

For  one-dimensional  problems,  Eqs  (7)  reduce  to 


Equation  (8)  is  the  grid  generation  equation  used  in  the  pres¬ 
ent  study.  It  is  solved  implicitly  by  an  optimized  Successive 
Over  Relaxation  (SOR)  technique.  One-sided  upwind  differences 
are  used  for  the  convective  term  based  on  the  results  of  a 
study  by  Ghia,  Hodge,  and  Hankey  (4).  These  authors  had  showed 
that,  when  the  control  function  became  large,  as  it  needs  to 
be  for  large  Reynolds  numbers,  and  a  central  difference  is 
used  for  the  convective  term,  the  solution  of  the  grid  equa¬ 
tion  becomes  nonmonotonic  and  an  oscillatory  solution  results. 
The  use  of  an  upwind  difference  was  shown  to  eliminate  this 
behavior . 

Adaptive  Grid  Criterion 

Any  finite-difference  representation  of  a  derivative  has 
truncation  error  associated  with  it.  This  truncation  error 
must  be  small  in  order  to  obtain  an  accurate  solution  of  the 
problem.  A  central-difference  representation  of  a  first  de¬ 
rivative  is 
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x 


f.  ,)/  2Ax  +  T.E. 
l-l 


where  the  truncation  error  (T.E.)  is  given  by 


T.E. 


f 

XXX 
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5.  xxxxx 


H.O.T. 


(9a) 


(9b) 


and  H.O.T.  stands  for  'higher  order  terms'.  Similarly,  a 


second-order  accurate  backward-difference  approximation  for 
a  first  derivative  is 


with 


f 


x 


T.E. 


=  (3f.  -  4fi_1  +  fi_2)/  2Ax  +  T.E. 


f 

XXX 


3 

+  f  +  H.O.T. 

4  xxxx 


(10a) 


(10b) 


In  the  past,  efforts  to  minimize  trucation  error  have 
consisted  of  reducing  the  grid  spacing  Ax.  This  does  indeed 
decrease  the  magnitude  of  the  truncation  error,  but  at  the 
expense  of  an  increased  number  of  grid  points  necessary  to 
cover  the  domain  and  therefore  an  increased  computational 
time.  For  many  problems,  this  may  not  be  insignificant,  es¬ 
pecially  for  higher  dimensions.  An  alternate  method  for  de¬ 
creasing  error  is  to  reduce  the  magnitude  of  the  higher  de¬ 
rivatives.  In  the  physical  plane,  one  has  no  control  over 
these  derivatives;  however,  in  the  computational  plane,  one 
could  gain  control  over  these  terms  if  the  transformation 
from  the  physical  plane  to  the  computational  plane  were  based 
on  the  reduction  of  these  derivatives. 

The  one-dimensional  transformation  of  a  first  derivative 
is  given  by  Eq  (2a) 

«x  •  £C  '  \  <2a) 

If  the  derivatives  in  the  computational  plane  are  expressed 


as  standard  central  differences  with  A£  taken  as  unity, 
then 


£x  •  <fi+l  -  fl-l>  '  2*5  +  T-E- 

(11a) 

with 

T.E.  -  -f?5?  /  6  x5  ♦  H.O.T. 

(lib) 

Thompson  (1:4-6)  gives  a  lengthy  argument  that  T.E.  must  be 
expressed  in  the  physical  plane  giving 


T.E. 


=  .  i  J£L  f 


.  i  x  f 
2  xx 


i  ; 
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f  + 

XXX 


H.O.T.  (He) 


Equation  (11c)  shows  the  T.E.  to  be  very  dependent  on  the 
grid  spacing.  In  short,  the  metrics  must  be  minimized  such 
that  their  expressions  in  Eq  (11c)  are  reduced  in  order  to 
reduce  T.E.  To  reduce  the  metric  error  completely,  however, 
would  be  to  eliminate  the  advantages  of  the  transformation,. 

On  the  other  hand,  if  one  considers  the  T.E.  of  f r  in 
the  computational  plane  (for  A?  =1)  then, 


T.E.  =  -f  /  6  +  H.O.T. 

S  '■a  S 

for  a  central-difference  representation  of  the  first  deriva¬ 
tive  and 


T.E.  *  -f  ^  /  3  +  H.O.T. 


(12b) 


for  a  second-order  backward  difference.  In  the  computational 
plane,  the  T.E.  depends  only  on  the  higher  order  derivatives. 
If  the  solution  in  the  computational  plane  were  a  second  de¬ 
gree  polynomial 

f.(C)  -  a2C2  +  ax5  +  aQ  (13) 

then  the  third,  fourth,  and  all  higher  order  derivatives  would 
be  identically  zero.  For  second-order-accurate  finite  dif¬ 
ferences,  the  truncation  error  would  be  eliminated,  no  matter 
what  the  grid  spacing,  and  an  optimum  grid  would  be  produced. 
The  number  of  grid  points  could  then  be  reduced  without  sacri¬ 
ficing  accuracy. 

In  general,  Eq  (13)  cannot  be  enforced  over  the  entire 
domain  because  there  will  inevitably  be  some  error  in  the 
solution  process  and,  if  Eq  (13)  is  not  satisfied  exactly, 
large  variations  in  the  grid  spacing  may  cause  significant 
errors  in  the  solution. 

It  may  only  be  necessary  to  enforce  Eq  (13)  locally. 

The  leading  term  of  the  local  truncation  error  is  proportion¬ 
al  to  the  third  derivative.  Let  it  be  approximated  by  a  cen¬ 
tral  difference,  giving 


T  E  -  f  -  f 


i+1 


(14a) 


i-1 


Again,  using  central  differences  to  evaluate  the  second  de¬ 
rivatives  gives 


13 


(14b) 


T.E.  . 
1 


'  i+2 


-  2f 


i  +  1 


+  2f 


i-1 


-  f i-2 


Therefore,  the  local  T.E.  at  a  point  is  approximately  depend¬ 
ent  only  on  its  immediate  neighbors.  In  this  context,  the 
entire  domain  can  be  divided  into  several  segments,  each  of 
which  satisfy  Eq  (13),  i.e. 


:i<5>  ' 


‘2l  5 


(15) 


At  each  grid  point,  a  second  order  least  squares  curve  fit 
is  performed  using  the  solution  at  the  previous  iteration  to 
calculate  the  constants.  The  number  of  data  points  NP  used 
in  the  calculation  may  be  specified  to  be  any  number  greater 
than  3  and  less  than  or  equal  to  the  total  number  of  grid 
points.  For  example,  if  NP  *  5,  then  the  constants  at  each 
grid  point  i  will  be  determined  by  the  previous  solution 
values  at  i±2,  i±l,  i  for  a  balanced  curve  fit.  At  the 
boundaries,  the  curve  fit  will  necessarily  be  unbalanced. 

The  grid  control  function  is  then  determined  as  described 
in  the  next  section.  If  NP  is  equal  to  the  total  number  of 
grid  points,  the  entire  domain  is  used  in  the  least  squares 
curve  fit,  producing  just  one  set  of  constants.  Because  the 
constants  are  based  on  the  solution  at  the  previous  iteration, 
the  grid  will  change  as  the  flow  solution  changes,  lagging  it 
by  one  iteration. 


14 


Grid  Control  Function  Evaluation 


The  previous  section  described  the  goal  of  the  adaptive 
grid.  In  this  section,  the  means  of  achieving  that  goal  is 
defined . 

Equation  (2a)  gives  the  transformation  of  the  first  de¬ 
rivative.  Solving  for  the  computational  derivative  gives 


f  , 


Differentiating  w.r.u.  £  gives 


(16) 


f 


+ 


f 


x 


(17) 


Using  the  grid  generation  equation,  Eq  (8),  and  Eq  (15),  an 
equation  of  the  form  F(P)  =  0  results  where 


F(P)  -  £xx  P£x  -  2  a2  /  x5  (18) 

for  which  a  Newton-Raphson  iteration  can  be  performed  to  solve 
for  the  root  P. 


pS*1  -  -  FS(P)  /  (|5)  S  (19> 

where  s  denotes  the  iterate  level.  Using  Eq  (18)  to  deter¬ 
mine  F(P)  and  —  gives 

a? 

Ps+1  =.  Ps  .  (f^  x^  -  Pfx  -  2a2/xc)/(-fx)  (20) 

At  this  point,  a  relaxation  factor  0  is  introduced  so  that 
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the  grid  movement  be  gradual  (9:321-322)  and  the  physical 
derivatives  are  replaced  by  their  computational  counterparts, 
Eqs  (2),  giving 


P 


s+1 


*  ps  +  e 


1 - 2a":£+  a 

L  C  2 


(21) 


where  0  <  9  <  1  and  f  in  the  last  term  has  been  replaced 

by  its  counterpart  from  Eq  (15)  for  consistency. 

All  values  at  s  are  known  for  the  previous  time  step 
or  iteration.  The  constants  a2  and  a^  are  evaluated  from 
the  least  squares  curve  fit  to  that  solution.  All  derivatives 
are  evaluated  as  second-order  central  differences  and  initially 
P  must  be  provided  as  input  to  start  the  solution.  Currently, 
no  formal  method  is  used  to  determine  9  so  it  must  be  deter¬ 
mined  by  experimentation. 

With  the  new  control  function  specified  a  new  grid  is 
then  determined  as  the  solution  of  Eq  (8).  For  comparison, 
the  grid  control  function  used  by  Freeman  (8)  is  given  by 
the  first  two  terms  on  the  right  hand  side  of  Eq  (21).  It  is 
referred  to  as  a  linear  method  because  it  results  from  speci¬ 
fying  the  solution  in  the  computational  plane  to  be  linear, 
therefore  producing  f  5  0.  It  is  felt  that  this  may  be 

SS 

over-constraining  the  problem  because  the  present  analysis 


shows  that  it  is  only  necessary  to  satisfy  the  condition 
f  -  constant,  (not  necessarily  0).  The  present  case  is 


referred  to  as  a  quadratic  method. 

One-Dimensional  Model  Problem 

The  viscous  Burgers'  equation  is  chosen  to  test  the  a- 
daptive  grid  method  because  it  is  typical  of  many  problems 
encountered  in  fluid  mechanics.  It  is  a  nonlinear,  second 
order  equation,  and  its  solution  produces  large  gradients  as 
the  Reynolds  number,  Re,  is  increased.  The  resulting  flow  can 
be  compared  to  a  boundary  layer  profile  of  thickness  propor¬ 
tional  to  ( 1/Re)  . 

The  nonconservative  form  of  Burgers'  equation  is 

Ut  +  UUx  -  Uxx  ' 
with  boundary  conditions 

U(-  °°  ,t)  =  1.0  U(0,t)  =  0.0 

and  initial  conditions 

U(x,0)  =  1.0  for  x  <  0  U(0,0)  =  0.0  (22c) 

The  analytical  steady  state  solution  is  given  by 

U(x)  *  -  tanh  (x  Re/2)  (23) 

Equation  (22a)  is  transformed  to  the  computational 


(22a) 


(22b) 


plane  by  use  of  Eqs  (2)  resulting  in 


V  + 

r 


U 


Re 


(24) 


An  optimized  Successive  Over  Relaxation  (SOR)  method  is 
used  to  solve  Eq  (24).  For  high  Reynolds  number  flows,  it 
has  been  shown  that  central  differencing  of  the  convective 
term  leads  to  an  oscillatory  solution  (4:23).  Therefore, 
second-order,  one-sided  upwind  differencing  is  used  for  the 
convective  term.  The  time  derivatives  are  expressed  as  a 
first-order  backward  differences  and  the  diffusion  term  is 
represented  as  a  second-order  central  difference.  The  metrics 
are  evaluated  as  central  differences. 


Solution  Procedure 

The  solution  algorithm  may  be  summarized  as  follows; 

The  computer  program  developed  to  solve  this  problem  using  the 
adaptive  grid  method  developed  here  is  listed  in  Appendix  B. 

1.  Provide  an  initial  guess  for  the  control  function  P. 

2.  Solve  Eq  (8)  for  the  grid  point  distribution. 

3.  Solve  Eq  (24)  to  obtain  the  flow  solution  for  the 
first  iteration. 

4.  Perform  a  least  squares  curve  fit  of  the  flow  solution 
to  determine  the  constants  and  a^. 

5.  Generate  a  new  grid  control  function  distribution  by 
solving  Eq  (21). 

6.  Repeat  Steps  2  through  5  until  steady  state  is 
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1 


reached . 


Steady  state  is  assumed  to  have  been  achieved  when  the  average 
difference  in  the  solution  between  two  successive  iterations 
is  less  than  a  specified  tolerance.  This  condition  is  given 
as 

\  r  |  U  n  -  u.n_1  I  <  £  (25) 

1  i=l  1  1 

where  I  is  the  total  number  of  grid  points,  n  denotes  the 
iteration  level,  and  t  is  the  tolerance. 


Error  Analysis 

For  this  study,  an  analytical  steady  state  solution  ex¬ 
ists,  thereby  providing  a  direct  means  to  determine  the 
truncation  error.  However,  in  general,  nc  such  solution  will 
be  available.  Another  measure  of  the  truncation  error  is 
defined  as  the  residual  error  R.E.  and  is  determined  as  the 


difference  between  the  computed  solution  and  the  solution  de¬ 
fined  by  Eq  (15).  The  local  error,  both  truncation  and  re¬ 
sidual,  is  defined  to  be  the  error  at  a  point.  The  global 
error  is  the  value  of  the  local  errors  averaged  over  the  en¬ 
tire  field ,  i  .  e  .  , 


R.E. 


U. 

l 


(26b) 


If  the  solution  is  such  that  the  R.E.  is  zero,  then  Eq  (15) 
would  be  satisfied  and  the  optimum  grid  would  be  defined. 
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III.  Discussion  of  Results 

The  adaptive  grid  method  developed  in  the  previous 
section  was  tested  by  solving  the  one-dimensional  Burgers' 
equation  given  in  Eq  (23).  The  solution  of  Burgers'  equa¬ 
tion  exhibits  an  increase  in  slope  at  the  right  boundary 
and  a  decrease  in  slope  in  the  left  boundary  as  Re  increas¬ 
es.  For  large  values  of  Re,  the  solution  of  Burgers'  equa¬ 
tion  is  comparable  to  boundary- layer  flows  except  that  the 
thickness  of  the  boundary  layer  is  proportional  to  (1/Re) 

i- 

instead  of  (1/Re)2.  An  essential  feature  of  any  mesh  is  a 
concentration  of  grid  points  in  regions  of  high  gradients 
in  the  flow  solution.  For  boundary-layer  flows,  a  good  rule 
of  thumb  is  to  have  5  to  10  grid  points  inside  the  boundary 
layer  thickness.  It  is  also  necessary  to  have  a  smooth  var¬ 
iation  in  the  distribution  of  grid  points  in  the  transition 
region  from  large  gradients  to  small  gradients  and  to  have 
enough  grid  points  in  all  regions  of  the  domain  so  as  to 
give  an  accurate  solution  of  the  problem.  These  features 
are  used  as  a  test  of  the  effectiveness  of  the  adaptive  grid 
method  developed  here. 

Basis  of  Results 

It  was  necessary  to  limit  the  value  of  the  control  func¬ 
tion  P  determined  by  Eq  (21)  because  early  results  produced 
very  large  magnitudes  of  P,  in  some  cases  resulting  in  a 
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double-valued  transformation  violating  the  maximum  principle. 
For  the  results  presented  here,  tha t  limit  was  set  at  2.0 
based  on  a  truncation  error  analysis  of  the  grid  generation 
equation,  Eq  (8).  For  a  uniform  distribution  of  P  =  2.0, 

I  *  21,  and  AC  *  1  analytical  solution  of  Eq  (8)  for  B.C.s 
x( 1)  =  -1  and  x(21)  =  0  gave  x(2)  -  -.0067.  The  correspond¬ 
ing  numerical  solution  using  a  second  order  upwind  difference 
for  the  convective  term  gave  x(2)  =  -.314.  Two  features  are 
evident  from  this  analysis.  First,  this  is  an  enormous  spac¬ 
ing  which  is  undesirable  for  accurate  solution  of  the  flow 
equations.  Second,  the  numerical  solution  is  very  different 
from  the  analytic  solution.  The  truncation  error  analysis 
showed  that  the  T.E.  is  proportional  to  exp(P).  The  limit 
value  of  2.0  was  chosen  because  the  solution  of  Eq  (22)  pro¬ 
duced  much  smaller  values  over  most  of  the  domain  and  local 
regions  of  P  =  2.0  could  be  handled.  An  additional  constraint 
was  imposed,  which  provided  a  more  uniform  distribution  of 
grid  points  near  the  wall  (x  =  0)  boundary.  This  constraint 
was  derived  from  the  transformed  Burgers'  equation,  Eq .  (24). 
Because  an  upwind  difference  was  used  for  the  convective 
term,  it  was  desirable  to  keep  the  value  of  the  coefficient 
multiplying  it  positive  adjacent  to  the  wall  in  order  to  use 
a  second-order  accurate  difference.  For  this  to  happen,  the 
value  of  P  had  to  be  less  than  the  product  U  •  Re  •  X  .  The 
only  place  where  this  expression  had  any  effect  was  near  the 
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wall  where  the  grid  spacing  was  very  small. 

For  this  study,  the  grid  was  defined  to  be  converged 
when  the  average  value  of  the  grid  speed  x  was  less  than 

T 

a  specified  minimum  e  .  The  final  converged  grid  depended 

xx 

on  several  factors: 

1.  The  initial  guess  of  the  control  function  P  which 
provides  the  initial  grid.  A  value  of  zero  was 
chosen  because  it  produces  a  cartesian  grid  of 
uniform  spacing.  In  this  way,  the  robustness  of 
the  grid  adaption  mechanism  could  be  demonstrated. 

2.  The  number  of  data  points  NP  used  in  the  least 
squares  curve  fit  which  determines  the  constants 
a 2  >  3^  and  3q . 

3.  The  grid  covergence  criteria  e  .  This  factor  was 

XT 

parameter  dependent  and  will  be  discussed  later  in 
this  section. 

4.  The  value  of  the  relaxation  factor  e  in  Eq  (8)  plays 
a  role  in  the  determination  of  e  .  Its  effect 

XT 

will  also  be  discussed  later  in  this  section. 

The  infinity  boundary  condition  for  the  velocity  U  was  speci- 
field  to  occur  at  x  5  -1.0.  For  small  values  of  the  Re,  U(-l) 
is  not  necessarily  equal  to  one,  therefore  it  was  determined 
from  the  analytic  steady  state  solution.  Results  are  pres¬ 
ented  here  for  21  grid  points  in  the  domain  -1.0  <  x  <  0.0 
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and  values  of  NP  =  5  and  NP  =21. 


Grid  Dependence  on  NP 

Figures  2  through  6  present  the  steady  state  solution 
and  the  converged  grid  for  NP  =21  and  for  Re  =  1 ,  10,  100, 
1000,  1500.  Large  slopes  in  the  grid  curve  indicate 
large  spacings  between  the  grid  points,  likewise  small  slopes 
in  the  grid  curve  indicate  small  spacings  between  the  grid 
points.  The  solution  in  the  transformed  plane  takes  on  a 
parabolic  shape  (as  expected)  for  the  smaller  values  of  Re. 

A  very  good  grid  is  produced,  placing  8  and  6  grid  points 
inside  the  boundary  layer  (B.L.)  for  Re  =  10  and  100,  re¬ 
spectively.  However,  as  Re  becomes  large,  the  solution  looks 
more  like  its  hyperbolic  tangent  form  in  the  physical  plane 
given  by  Eq  (23).  This  is  because  the  velocity  U  is  equal 
to  1.0  for  all  but  very,  very  small  values  of  x,  because 
the  boundaries  are  fixed,  and  because  a  quadratic  function 
is  a  very  poor  fit  to  the  hyperbolic  tangent  function. 

Also  note  the  step  in  the  grid  point  distribution  for 
Re  =  1000  and  1500.  This  step  is  due  to  the  last  term  of 
Eq  (21).  The  coefficients  &2  an^  ai  were  generally  of  op¬ 
posite  signs  and  a^  was  generally  an  order  of  magnitude 
larger  than  a2«  F or  small  values  of  £  the  sign  of  this  term 
was  determined  by  the  sign  of  a^;  however,  as  C  increased 
in  value  to  approximately  8  or  9,. the  2a9?  term  dominated. 
Also,  for  small  values  of  5  »  this  entire  term  dominated 
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the  first  term,  Ur-/Ur,  which  for  large  Re  was  equal  to  zero 

S  *9  ~3 

due  to  (J  being  constant  in  value.  For  larger  values  of 
the  U-r/Ur  term  dominated.  The  net  effect  was  to  produce 
the  grid  control  function  distribution  shown  in  Figure  7 
for  Re  *  1500.  Positive  values  of  P  cause  the  grid  points 
to  move  to  the  right,  X  =  0  boundary,  whereas  negative 
values  of  P  cause  movement  in  the  opposite  direction.  Where 
P  changes  sign,  the  grid  points  on  either  side  are  forced 
together,  therefore,  producing  the  step  in  Figures  5  and  6. 
Because  of  this  step,  several  grid  points  were  placed  in 
a  "pseudo  gradient"  leaving  only  a  few  grid  points  for  the 
real  gradient.  Only  two  grid  points  were  inside  the  bound¬ 
ary-layer  thickness  of  (1/Re)  for  Re  *  1000  and  for  Re  =  1500, 
only  1  grid  point  was  inside  the  boundary  layer. 

This  effect  is  eliminated  when  the  least  squares  curve 

fit  is  performed  only  locally,  i.e.,  for  small  values  of  NP. 

Figure  8  through  12  present  results  for  NP  =  5  for  the  same 

Reynolds  numbers  as  before.  For  Re  =  1 ,  the  grid  did  not 

move  significantly  from  its  original  position.  The  pseudo 

gradient  is  no  longer  present  for  the  larger  values  of  Re 

because  it  is  much  easier  to  fit  a  quadratic  function  for 

only  a  few  local  data  points  than  for  the  entire  domain.  The 

number  of  grid  points  inside  the  boundary  layer  is  7  for  Re 

*  100,  4  for  Re  =  1000,  and  6  for  Re  =  1500.  When  c  was 

x 

T 

reduced  to  .001,  7  points  were  placed  inside  the  B.L.  for 
Re*1000.  Figure  13  presents  a  grid  control  function  dis¬ 
tribution  for  these  cases  (NP*5)  showing  a  much  better  dis- 


tribution  than  the  NP=21. 


The  number  of  data  points  used  in  the  least  squares 

curve  fit  used  to  determine  the  constants  a2>  a^,  and  ag 

has  a  significant  effect  on  the  final  converged  grid.  Figure 

14  demonstrates  as  NP  is  increased  the  converged  grid  looks 

more  and  more  like  that  for  NP  =21.  Even  if  the  value  of  the 

grid  convergence  criteria  e  is  decreased,  the  final  grid  will 

xt 

exhibit  this  step  for  the  larger  values  of  NP. 

Grid  Dependence  on  e„ 

It  is  evident  that  the  value  of  the  grid  convergence 
criteria  is  a  very  important  factor  in  the  shap  of  the  final 
grid  and  in  the  accuracy  of  the  flow  solution.  As  the  para¬ 
meters  (e.g.,  Re,  NP,  etc.)  are  changed,  it  becomes  difficult 
to  predict  the  value  of  e  which  will  produce  the  best  grid 

At 

due  to  the  highly  oscillatory  nature  of  the  grid  speed  vari¬ 
ation  with  time.  Figure  15  shows  a  representative  time  his¬ 
tory  of  the  average  grid  speed  demonstrating  this  feature. 

If  too  small  a  value  of  e  x  was  chosen,  the  grid  did  not  con- 

T 

verge  in  a  reasonable  number  of  iterations.  For  comparison, 
the  solution  of  this  model  problem  on  a  static  grid  required 
between  5  to  17  iterations  to  converge  to  steady  state. 

Steady  state  solution  on  an  adaptive  grid  required  between  15 
to  30  iterations,  on  the  average.  However,  if  e x  was  too 

T 

small,  grid  convergence  was  not  attained  even  after  50  itera¬ 
tions.  If  too  large  a  value  of  ex  was  chosen,  the  converged 
grid  did  not  fully  resolve  the  gradients  in  the  flow  solution. 
Figure  16  demonstrates  this  behavior. 


Figures  17  through  19  present  a  steady  state  error  analy¬ 


sis  based  on  £x  for  two  different  values  of  Re  and  two  values 

T 

of  NP.  The  key  quantity  in  these  figures  is  the  maximum  trun¬ 
cation*  error  ,  T.E.^.  It  might  be  expected  that,  if  the  value 

of  e  was  smaller,  the  solution  would  be  improved.  Figure 
x  T 

18  and  19  show  this  to  be  true.  However,  Figure  17  shows  the 

contrary.  At  this  point,  no  definite  value  can  be  proposed 

for  an  arbitrary  set  of  conditions.  The  key  quantity  in 

these  figures  is  the  maximum  truncation  error,  T.E.M.  Ot 

might  be  expected  that  if  the  value  of  e  was  smaller,  the 

x 

solution  would  be  improved. 


Grid  Dependence  on  9 

The  value  of  the  relaxation  factor  9  has  an  effect  on 
on  the  final  grid  point  distribution,  although  not  as  much 
as  the  other  parameters.  As  0  was  increased,  more  of  the 
grid  movement  generated  by  Eq  (21)  was  allowed  to  be  accom¬ 
plished  at  each  iteration,  i.e.,  the  grid  speed  increased. 
This  may  be  good,  especiallyf or  larger  Re  where  more  grid 
movement  is  necessary  to  resolve  the  flow  gradients.  Along 
with  this  increased  grid  speed  comes  an  increase  in  the  amp¬ 
litude  of  the  oscillation  in  the  grid  speed  convergence,  dis¬ 
cussed  previously.  The  smaller  the  value  of  e,  the  more  uni¬ 
form  the  grid  movement  became.  Figure  20  shows  a  representa¬ 
tive  effect  of  9  on  the  steady  state  error.  In  general,  a 
better  solution  was  achieved  as  9  was  decreased,  however, 


there  is  a  drawback.  In  general,  the  smaller  -a  was,  the 
longer  convergence  took. 


Comparison  with  Other  Methods 

The  effectiveness  of  the  adaptive  grid  method  developed 
in  this  thesis  is  now  compared  to  the  adaptive  grid  method 
of  Freeman  (8)  and  to  a  static  grid.  An  attempt  was  made 
to  compare  the  different  methods  on  the  same  basis;  a  key 
feature  of  the  comparison  is  the  resolution  of  the  flow  grad¬ 
ients.  A  measure  of  this  is  the  number  of  grid  points  in  the 
boundary  layer  (1/Re).  Therefore,  for  high  Re,  the  static 
grid  was  exponentially  stretched  in  order  to  place  a  similar 
number  of  grid  points  inside  the  boundary  layer  (1/Re).  As 
many  parameters  as  possible  wer  kept  the  same,  however,  in 
many  instances,  this  was  not  possible. 

Table  I  summarizes  the  comparison  of  the  three  methods 
for  Reynolds  numbers  of  10,  100,  1000,  and  1500.  The  key 

parameter  to  note  is  the  maximum  truncation  error  (T.E.  ). 

M 

Freeman's  method  is  referred  to  as  the  linear  method  and  the 
present  method  is  referred  to  as  quadratic.  For  NP  =  5, 
the  adaptive  methods  are  very  similar,  which  is  demonstrated 
in  Figure  21.  The  values  T.  E.  for  the  adaptive  grids  com¬ 
pare  very  well  with  the  T.E.  static  of  the  grid.  Note,  how¬ 
ever,  that  the  number  of  grid  points  is  substantially  increased 
for  the  static  grid  for  Re  =  10  and  100  in  order  to  resolve 
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TABLE  I 


Error  Comparison 
6  =  0.25  1=21  NP  =  5 


T.E.  max 

T.E 

:.  ave 

Re 

Fixed 

Adaptive  Grid 

Fixed 

Adaptive  Grid 

Grid 

Linear 

Quadratic 

Grid 

Linear 

Quadratic 

10 

.0021° 

.0093 

.0111 

.0007 

.0055 

.0042 

100 

. 0 113b 

.0198 

.0115 

.0005 

.0103 

.0053 

1000 

. 0169c  ! 

t 

.0227 

.0291 

.0049 

.0084 

.0087 

1500 

. 0160  c 

.0248 

.0365 

.0046 

.0078 

.0073 

a  1  =  51 

b  I  =  201 

c  I  =  21  P  =  0.5 


f 
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the  flow  gradient.  For  large  Re,  the  static  grid  was  expon¬ 
entially  stretched  to  resolve  the  gradient.  This  amounts  to 


a  constant  value  of  P  in  the  grid  generation  equation,  Eq  (8). 
This  stretching  was  possible  for  this  model  problem  because 
the  characteristics  of  the  solution  were  known,  however,  in 
general,  these  features  will  not  be  known  a  priori  and  a 
stretching  of  this  sort  would  not  be  possible.  One  would 
then  have  to  resort  to  a  large  number  of  grid  points  in  order 
to  resolve  all  the  features  of  the  flow. 

Another  important  factor  concerning  the  viability  of 
any  new  develpment  is  the  relative  computation  time  required 
to  solve  the  problem.  This  study  was  performed  on  a  CDC 
Cyber  175  computer.  For  the  static  grid  the  solution  of 
Burgers'  equation  required  an  average  of  2.03  x  10  ^  cp 
seconds  per  grid  point  per  iteration.  Freeman's  linear  adap- 

_3 

tive  grid  method  required  an  average  of  3.88  x  10  cp 

_3 

seconds  and  the  present  quadratic  method  required  4.31  x  10 

-3 

seconds  for  NP  =*  21  and  4.67  x  10  cp  seconds  for  NP  =  5. 

The  present  method  required  between  2  and  5  cp  seconds  to  r 
reach  steady  state  for  21  grid  points  compared  to  between  .2 
and  .5  cp  seconds  for  21  grid  points  for  the  static  grid. 
However,  the  static  grid  with  201  grid  points  required  3.066 
cp  seconds. 
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IV.  Conclusions 


The  results  of  this  study  demonstrate  the  benefits  of 
using  an  adaptive  grid  in  the  solution  of  fluid  dynamics 
problems.  There  were  two  primary  goals  of  the  adaptive  grid 
method  develped  in  this  thesis:  first,  to  provide  adequate 
resolution  of  the  high  gradient  regions  in  flow;  and  second, 
to  produce  an  optimum  grid  such  that  the  truncation  error 
would  be  eliminated.  The  first  objective  was  met.  The 
method  does  a  very  good  job  of  concentrating  grid  points  in 
the  physical  plane  in  high  gradient  regions,  where  large 
truncation  errors  generally  occur.  It  showed  robustness  in 
that  the  gradients  were  resolved  given  an  initially  constant 
spaced  cartesian  grid  and  without  any  a  priori  knowledge  of 
the  flow.  However,  it  is  not  as  robust  as  hoped,  because 
it  depends  on  the  input  of  several  parameters  which  presently, 
can  only  be  determined  by  experimentation.  In  general,  the 
following  statements  can  be  offered  for  the  determination 
of  these  parameters  based  on  the  results  of  this  study. 

1.  The  number  of  data  points  used  in  the  least  squares 
curve  fit  should  be  decreased  as  Re  increases. 

Using  the  entire  field  for  the  curve  fit  was  suc¬ 
cessful  only  for  smaller  values  of  Re.  As  the  num¬ 
ber  of  points  is  decreased,  the  method  approaches 
that  of  the  linear  method  used  by  Freeman  (8). 
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2.  Large  values  of  the  relaxation  factor  6  result 

in  greater  grid  speeds  and  less  convergence  times. 
Smaller  values  result  in  decreased  grid  speeds  and 
larger  convergence  times.  However,  the  smaller 
values  of  6  produce  better  steady  state  results  due 
to  the  dampening  of  the  oscillations  that  occur  in 
the  grid  speed. 

3.  At  this  point,  no  conclusion  can  be  drawn  about  the 
grid  convergence  criteria  e  .  It  appears  to  be 

'  T 

very  problem  dependent  and  can  only  be  determined 
by  experimentation. 

4.  The  computation  time  required  for  solution  on  an 
adaptive  grid  is  greatly  increased  over  that  on  a 
static  grid.  However,  without  a  priori  knowledge 
of  the  flow  solution,  a  great  number  of  grid  points 
is  required  for  the  static  grid,  putting  the  two 
methods  on  equal  ground  for  this  one-dimensional 
problem . 

The  second  goal  of  the  thesis  was  achieved  only  partially. 
The  method  produces  satisfactory  results;  however,  it  did  not 
produce  a  truly  optimum  grid.  The  truncation  errors  are  still 
on  the  same  order  as  those  produced  by  a  static  grid  with  e- 
qual  grid  spacing  or  with  coordinate  stretching. 
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V.  Recommendations 


The  results  of  this  study  indicate  the  advantages  of 
the  use  of  a  solution  adaptive  grid  for  the  numerical  solu¬ 
tion  of  partial  differential  equations.  Further  investi¬ 
gation  of  the  present  method  is  required  to  determine  the 
reason  the  truncation  errors  were  not  minimized,  as  ex¬ 
pected. 

One  possible  source  of  error  resided  in  the  Newton 
iteration  technique  used  to  determine  the  grid  control 
functions.  It  is  suggested  that  a  fixed  Newton  iteration  be 
used;  the  denominator  should  consist  of  only  one  term,  not 
two . 

A  second  source  of  error  is  associated  with  the  grid 
spacing.  Large  errors  in  the  solution  tended  to  occur  where 
large  grid  spacings  occurred,  especially  if  the  spacing  was 
three  or  more  times  the  value  of  the  smaller  grid  point. 

The  present  method  may  need  to  be  modified  to  restrict 
the  size  of  the  grid  spacing. 


► 
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Re  -  1 


0-0.25  «  -.0001  NP  -  21 

xr 


Fig  2  Velocity  Profile  &  Grid  Re  =  1  NP  =  21 
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Velocity,  U  ;  Physical  Coordinate, 


Fig  3  Velocity  Profile  &  Grid  Re  =  10  ,  NP  «  21 
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Velocity,  U  ;  Physical  Coordinate, 


sical  Coordinate 


Velocity,  U  ;  Physical  Coordinate 
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Computational  Coordinate,  £ 


;-'ig  6  Velocity  Profile  &  Grid  Re  =  1500  ,  NP  =  21 
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A 


Grid  Control  Function 


Velocity,  U  ;  Physical  Coordinate, 


Fig  8  Velocity  Profile  &  Grid  Re  ■  1  ,  NP  *  5 
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Velocity,  U  ;  Physical  Coordinate 


Re  -  10 


..001 


NP 


e-  1.0 


Fig  9  Velocity  Profile  &  Grid  Re  ■  10  ,  NP  » 
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Velocity,  U  ;  Physical  Coordinate, 


Fig  10  Velocity  Profile  &  Grid  Re  *  100  ,  NP  *  5 
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Velocity,  U  ;  Physical  Coordinate, 


Fig  11  Velocity  Profile  &  Grid  Re  =  1000  ,  NP  -  5 
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Fig  12  Velocity  Profile  S>  Grid  Re  *  1500  ,  NP  *  5 
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Re  =  1500 


9  ~  1.0  €  =.005  NP  =  5 
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i  i  ■  i  i 

0  5  10  15  20 

Computational  Coordinate,  £ 

Fig  13  Grid  Control  Function  Distribution,  Re  =  1500,  NP  =  5 
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Physical  Coordinate, 


Re  =  1000 


$=  0.25 
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Average  Or  id  Speed, 


Time,  t  j 

Fig  15  Grid  Speed  Time  History  I 
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Fig  18 


Fig  19  Solution  Error  Depender 


Velocity,  U  :  Physical  Coordinate, 


Fig  21  Grid  Comparison  with  Linear  Mothod 
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Agoendi_x_B 


PROGRAM  LISTING 

PROGRAM  BURG 1 5C (I NPUT , OUTPUT , TAPE6 , TAPE7 , TAPES ) 

cmmmmmmmmmmmmmmmmmmmmm 

C  UNSTEADY  SOLUTION  OF  THE  1-D  BURGER’S  EQUATION  BY  AN  * 

C  OPTIMIZED  SOR  MITHOD  COUPLED  WITH  A  GRID  OPTIMIZATION  * 

C  ROUTINE  BASED  ON  A  TRUNCATION  ERROR  ANALYSIS.  THE  * 

C  PROCEDURE  IS  NOT  SELF-STARTING.  THUS  REQUIRING  * 

C  AN  INITIAL  GUESS  FOR  THE  GRID  GENERATION  CONTROL  * 

C  PARAMETER,  P.  * 

cmmmmmmmmmmmmmmmmmmmmm 

COMMON  /A/ IMAX, U (51) , ZI (51 ) , X (51) ,P(51> , DX (51 ) , RE 

COMMON  /B/KGRID, GRIDACC, XMIN, XMAX 

REAL  DXX (51) , XT (51) , XN(51) ,UN(51) ,E(51 ) , G (51) , C(3) 

REAL  LBC, RBC, XMIN, XMAX, ALPHA, BETA, GAMMA, DX2, A, B, USTAR, W, 

#  UOLD, SUM, ERRMAX , ERRAVE, DIFF, GRIDACC. SOLNACC, DU, DUU, DELT, T 

#  XTAVE, XTACC.CPI , CPF, CPU, PI, ZZ (1 1) , UU( 11 ) 

INTEGER  I , K, N, NT, KGRID, KSQLN 

PI«3. 141592654 

ccmmm  read  input  data  mmmmmmmmmmmm 

READ (7,  *)  IMAX, KGRID, KSOLN, LBC, RBC, XMIN, XMAX, NP 
READ ( 7 , * )  RE , GR I DACC , SOLNACC , DELT , NT , THETA , XT ACC 
READ (7, t) <P(I) , 1*1, IMAX) 

CALL  DATE(ADATE) 

CALL  TIME(ATIME) 

WRITE (6, 45)  ADATE, ATIME 

WRITE (6, 50)  IMAX, NT, RE 

WRITE (6, 51)  LBC, XMIN, RBC, XMAX 

WRITE (6, S3)  KSOLN, SOLNACC, KGRID, GRIDACC 

WRITE (6, 54)  DELT, XTACC, THETA, NP 

WRITE (8, 200)  ADATE, ATIME 

WRITE (8, 205)  IMAX, NT, RE 

WRITE(8,210)  LBC, XMIN, RBC, XMAX 

WRITE (8, 215)  KSOLN, SOLNACC, KGRID, GRIDACC 

WRITE (8, 210)  DELT, XTACC, THETA, FLOAT (NP) 

ccmmm  set  initial  conditions  mmmmmmmm**m 

T-0.0 

N-0 

DO  10  1*1, IMAX 
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ZI <I) "FLOAT (I) 

X (I)*XMIN+(I-1) t (XMAX-XMIN) / ( IMAX-1 ) 

Ud)*1.0 
10  CONTINUE 

CALL  USTART 

CALL  UPSET (’POLYNOMIAL*, 2.0) 

XTAVE-1.0 

CCtttltUlt  SET  BOUNDARY  CONDITIONS  tttttttttUtntttUUtttttU 

U(l>*-TANH(RE*XU)t0.3> 

U( IMAX) — TANH (REIX  < IMAX) *0. 5) 

XT(1)*0.0 
XT (IMAX) -0.0 

CCUtttttt  CALCULATE  'OLUTION  FOR  EACH  TIME  STEP.  N  ************ 

CALL  SECOND (CPI. 

990  N-N+l 

T-T+DELT 

IF  (N  .ST.  NT)  GO  TO  1000 
K*0 

DO  103  1*1, IMAX 
XN(I)-X(I) 

103  UN(I)-Ud) 

WRITE <6, 63)  T 
WRITE(B,210)  T 

cctmtm  calculate  grid  ttmmmmmmmmmmmm 

IF  (N  .LE.  3  )  GO  TO  24 
IF  (ABS(XTAVE)  .LE.  XTACC)  GO  TO  23 

24  IUPWND-0 

CALL  GRIDSOR ( IUPWND) 

25  CONTINUE 

ccmmtt  write  grid  data  mm*mm*mm***mmm*«m 

XTAVE-0.0 
WRITE (6, 52) 

DO  20  1*2, IMAX-1 
DXd)*0.5*(X<I+l)-X(I-i) ) 

DXXd)*Xd+l)-2.0*X(I)+X(I-l) 

XT  d ) ■ ( X ( I ) -XN  d ) ) /DELT 
XTAVE*XTAVE+XT (I) 

WRITE (6, 38)  I,  Xd)  ,DXd) ,  DXXd) ,  XT(I)  ,P(I) 

20  CONTINUE 

XTAVE-XTAVE/dMAX-2) 

WRITE(6,59)XTAVE 

IF  (N  . GE.  NT-3)  XTAVE-0.0 

ccmtm*  k  is  the  sqr  iteration  loop  index  mmmmmmt 
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DO  900  K-l.KSQLN 
ERRMAX=0. 0 
SUM«0.0 

cc*****m  i  is  the  space  loop  index  *************************** 

DO  110  I«2, IMAX-1 
DX2-DX(I)*DX(I) 

ALPHA* 1 . 0/ (REXDX2) 

BETA- CU < I > -XT ( I )  ♦  DXX ( 1> (ALPHA) /DX ( I) 

C  USE  A  FIRST  ORDER  DIFFERENCE  NEXT  TO  THE  BOURNDARY  * 

C  POINTS  AND  A  SECOND  ORDER  DIFFERENCE  INTERIOR  * 

C  TO  THESE  POINTS  FOR  THE  UPWIND  CONVECTIVE  TERM  * 

C  THE  UPWIND  VALUE  IS  BASED  ON  THE  VALUE  OF  BETA  I 

c*************************************************************** 

IF  (BETA  .GE.  0)  THEN 

IF  (I  .EQ.  2  .OR.  I  .EQ.  ( IMAX-1 >>  THEN 
DU— U(I-l) 

GAMMA  ■  BETA  +  2.0* ALPHA  +  1.0/DELT 
A  «  ALPHA/GAMMA 
B  ■  (ALPHA  ♦  BETA) /GAMMA 
ELSE 

DU  -  0.5*U(I-2)  -  2.0*U(I-1) 

GAMMA  -  1.5IBETA  +  2.0*ALPHA  ♦  1.0/DELT 
A  •  ALPHA/GAMMA 
B  -  (ALPHA  +  BETA) /GAMMA 
END  IF 
ELSE 

IF  (I  .EQ.  2  .OR.  I  .EQ.  (IMAX-1))  THEN 
DU  ■  U(I+1) 

GAMMA  ■  -BETA  +  2.0*ALPHA  +  1.0/DELT 
A  *  (ALPHA  -  BETA) /GAMMA 
B  *  ALPHA/GAMMA 
ELSE 

DU  -  -0. 5*U ( 1+2)  ♦  2. 0$U( I+i ) 

GAMMA  -  -1 . 5*BETA  +  2.0*ALPHA  +  1.0/DELT 
A  -  (ALPHA  -  BETA) /GAMMA 
B  -  ALPHA /GAMMA 
END  IF 
END  IF 

CCS**!****  CALCULATE  THE  GAUSS-SEIDEL  VALUE  *X****«*I**I**$*«|** 

USTAR* (UN ( I ) /DELT+ALPHA* (U  < 1+ 1 >  +U  < 1-1 > ) -BETAtDU) /GAMMA 

CC*****t*t  CALCULATE  THE  OPTIMUM  RELAXATION  FACTOR  ************* 

CALL  NOPT(A,B,IMAX,W> 

IF (N  .EQ.  1  .AND.  K  .EQ.  1)  W*1.0 
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C*m  CALCULATE  SOLUTION  &  ERROR  BETWEEN  SOR  ITERATIONS  ******** 
UQLD  ■  U(I) 

U<I>  ■  U(I)  +  W* (USTAR  -  U(I) ) 

DIFF  «  ABS(U(I)  -  UOLD) 

SUM  -  SUM  +  DIFF 
110  CONTINUE 

CCM******  CHECK  FOR  ITERATION  CONVERGENCE  ********************* 
ERRAVE  -  SUM/IMAX 

CC  IF  (ERRAVE  .LE.  SOLNACC)  GO  TO  910 

900  CONTINUE 

910  WRITE (6, 79) 

«***»  CALCULATE  NEW  GRID  GENERATION  CONTROL  FUNCTION  P  ******** 

IF  (N  .LE.  3)  GO  TO  34 
IF  (ABS(XTAVE)  .LE.  XTACC)  GO  TO  35 

34  CALL  NEWP ( THET A , C , NP , G ) 

35  CONTINUE 

CC**t****«  CALCULATE  ANALYTIC  STEADY  STATE  SOLUTION  ************ 

WRITE <6. 87) 

ERR1-0.0 
SUM1-0.0 
ERR2-0.0 
SUM2-0.0 
ERR3-0. 0 
SUM3-0.0 
DO  120  1-1, IMAX 
E(I)— TANH<RE<X<D  *0.5) 

DIF1»A8S<U(I)-E<I)) 

SUM1-8UM1+DIF1 
DIF2-ABS (U ( I ) -6 ( I ) ) 

SUM2-SUM2+DIF2 
DIF3-ABS (G  < I ) -E  < I ) ) 

3UM3-SUM3+DIF3 

WRITE (6, 58) I,X(I),U(I),E(I) , DIF1, G(I) 

WRITE (8, 223) X (I) ,XT(I),P<I)»U<I),G<I),E(I) 

120  CONTINUE 

ERR A V l "SUM 1 / I M A X 
ERRAV2-8UM2/IMAX 
ERRAV3-SUM3/IMAX 

WR I TE  <  6 , 80 ) K- 1 , ERRAV 1 , ERR A V2 , ERR AV3 
WRITE (8, 230) K-l , ERRAV 1 , ERRAV2, ERRAV3 

CM***  CALCULATE  AVERAGE  ERROR  BETWEEN  SOR  ITERATIONS  ********** 

SUM-0.0 
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DO  180  1=1 , IMAX 

DIFF=ABS(U(I)-UN(I>) 

SUM-SUM+DIFF 
180  CONTINUE 

ERRAVE*SUM/ IMAX 
WRITE (6, 82) ERRAVE 
WRITE (8, 210) ERRAVE 
IF  (ERRAVE  .GE.  SOLNACC)  GO  TO  990 

CALL  SECOND (CPF) 

CPU-CPF-CPI 
WRITE (6, 95) CPU 
WRITE (8, 235) CPU 
CALL  UEND 
STOP 

200  FORMAT (T5,  A10, 5X, A10) 

205  FORMAT  <T5,2I5,E12.4) 

210  FORMAT (T5,4E12. 4) 

215  FORMAT (TS,  IS,  E12. 4, 15,  £12. 4) 

225  FORMAT (6E11.4) 

230  FORMAT (T5, IS. 3E12. 4) 

235  FORMAT (T3,E12. 4) 

45  FORMAT  <’ 1’ , /,  T5,  A10, 2X, A10) 

50  FORMAT (/, T10, ’ 1-D  BURGERS  EQ.  SOLN  USING  OPTIMIZED  SOR  METHOD' 
*,’  L  OPTIMUM  GRID’ , //, T25, ’********** INPUT  DATA**********’ , /3X, 

OF  GRID  PTS  »’ , 13, 3X, ’#  OF  TIME  STEPS  ’,I3,3X,’RE  #  *’,E9.3> 

51  FORMAT  (3X, ’U  ■’,F5.2,’  9  X  »’,F5.2,3X, ’U  »’,F5.2,’  9  X  «’,F5.2> 

53  FORMAT (3X, ’ 9  OF  SOLN  ITER.  , 13, 3X, ’ SOLN  ACC  CK 
#E9.3,/,3X,’#  OF  GRID  ITER.  »’ , 13, 3X, ’GRID  ACC  CK  »’,E9.3> 

54  FORMAT (3X,* DELTA  TIME  ■’ , F5. 2, 3X, ’ GRID  SPEED  CK  «’,E9.3,/, 

M3X , ' NEW  P  RELAXATION  FACTOR: ’ ,F4. 2, 3X, ’#  OF  PTS  IN  ’, 

#’ CURVE  FIT:’, 13,/) 

52  FORMAT (/,T25, 'GRID  CALCULATION' ,/, T4, ' ZI ', T14, ' X’ , T25, 

9’DX' ,T38, ’DXX’ ,T50, 'XT' ,  T62, ’P’ ,/) 

58  FORMAT (15, 8 (F12. 7)) 

59  FORMAT (/,T5,’ AVERAGE  VALUE  OF  GRID  SPEED  «’,F12.7> 
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65  FORMAT <//,T30. ’TIME  =’,F6.3./> 


79  FORMAT (//) 

30  FORMAT <T5. 'SOR  ITERATION  , I5.5X, ’ AVERAGE  ERROR  ANALYSIS’./. 
#T30, ’ COMPUTED  &  ANALYTICAL  « ’ , F13.8, /, T30, ’ COMPUTED  &  FIT  t’, 
#8X, F12. 8, /,T30, ’ANALYTICAL  &  FIT  i’,3X,F13.8> 

82  FORMAT <T5,’AVE  ERROR  BETWEEN  TIME  STEPS  »’,F12.B> 

87  FORMAT (/, T30, ******  ANALYTIC  SOLUTION  *****’,// 

#T4,’ZI*,T14,*X*, T25, * U’ , T33, *U  EXACT* , T47, ’ DIFF' , T60, ’FIT’ ,/> 

90  FORMAT*/, T20, ’LEAST  SQUARES  FIT  TO  SLUTION’ , /, T5, ’ U(ZI )  »  ’, 
#F10. 6, ’  *  ZI**2  +  * ,F10.6, *  *  ZI  +  ’,F10.6,/> 

95  FORMAT (/,T5, ’TOTAL  CPU  TIME  :’.  E17.7) 

98  FORMAT </, TS, ’MAX  ITERATIONS  EXCEEDED !!*!!',/) 

1000  WRITE (6, 98) 

CALL  SECOND (CPF) 

CPU*CPF-CPI 
WRITE <6, 95) CPU 
WRITE <8, 235) CPU 

STOP 

END 


SUBROUTINE  WOPT(A,B, IMAX,W) 

PI-3. 141592654 

PJ»/’.0*SQRT(A*B)*0.99 

W=2.0/ (1.0+SQRT  (1.0-PJ**2)  > 

RETURN 

END 


SUBROUTINE  GRIDSOR* IUPWND  ) 

c*t****tt*tt**t**mmtt****mtmt****t*t*m***m****m** 

C  THIS  SUBROUTINE  CALCULATES  THE  SOLUTION  OF  THE  1-D 
C  GRID  GENERATION  EQUATION  BY  THE  THOMAS  ALGORITHM. 

C  IF  IUPWND  «  1,  THEN  IT  USES  2ND  ORDER  UPWIND  DIFFERENCES 
C  FOR  THE  CONVECTIVE  TERM.  IF  IUPWND  ■  0,  IT  USES  A  2ND 
C  ORDER  CENTRAL  DIFFERENCE  FOR  THE  CONVECTIVE  TERM 
c************************************************************** 
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COMMON  /A/IMAX. U (51 ) ,  ZI (51) . X (51) , P  (31) , DUMMY (51 ) . RE 
COMMON  /B/KGRID. GRIDACC. XMIN. XMAX 
WRITE (6. 36) 

X  < 1 ) -XMIN 
X ( IMAX) »XMAX 
K  ■  0 

CCimm*  K  IS  THE  SOR  ITERATION  LOOP  INDEX  ******************* 

999  K  «  K+l 

IF  (K  .GE.  KGRID)  GO  TO  200 

ERRMAX-0.0 

SUM-0.0 

cc ********  i  is  the  space  loop  index  *************************** 

DO  110  1=2. IMAX-1 
ALPHA- 1.0 
BETA— P  <  I ) 

cttttttxttttttttttttttttttttttttttittttttttxtttttttinttxtttxtttt 


USE  a  FIRST  ORDER  DIFFERENCE  NEXT  TO  THE  BOURNDARY  * 
POINTS  AND  A  SECOND  ORDER  DIFFERENCE  INTERIOR  * 
TO  THESE  POINTS  FOR  THE  UPWIND  CONVECTIVE  TERM  * 
THE  UPWIND  VALUE  IS  BASED  ON  THE  VALUE  OF  BETA  * 


c*************************************************************** 

IF  (BETA  .GE.  0)  THEN 
IF  (I  .EQ.  2  .OR.  I  .EQ.  (IMAX-1) >  THEN 
DX— X(I-l) 

GAMMA  -  BETA  +  2. 01 ALPHA 
A  -  ALPHA/GAMMA 
B  -  (ALPHA  +  BETA) X GAMMA 
ELSE 

DX  -  0. 5*X < 1-2)  -  2.0*X(I-1) 

GAMMA  ■  1.5IBETA  +  2. Of ALPHA 
A  -  ALPHA /GAMMA 
B  =  (ALPHA  +  BETA) /GAMMA 
END  IF 
ELSE 

IF  (I  .EQ.  2  .OR.  I  .EQ.  (IMAX-1))  THEN 
DX  ■  X(I+1) 

GAMMA  -  -BETA  ♦  2.0* ALPHA 
A  -  (ALPHA  -  BETA) /GAMMA 
B  -  ALPHA/GAMMA 
ELSE 

DX  -  -0.5*X(I+2)  +  2.0<X(I+1) 

GAMMA  -  -1,3*BETA  ♦  2.0*ALPHA 
A  -  (ALPHA  -  BETA) /GAMMA 
B  -  ALPHA/GAMMA 
END  IF 
END  IF 
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cc ********  CALCULATE  THE  GAUSS-SEIDEL  VALUE  ******************** 


XSTAR  =  (ALPHA* (X ( 1+1 ) +X < 1-1 > ) -BETA*DX) /GAMMA 

CC********  CALCULATE  THE  OPTIMUM  RELAXATION  FACTOR  ************* 

CALL  W0PT(A,8,IMAX,W> 

IF  <K  ,EQ.  1)  W-1,0 

C*»***  CALCULATE  SOLUTION  &  ERROR  BETWEEN  SOR  ITERATIONS  ******* 
XOLD  »  X(I) 

X(I)  -  X(I)  +  W* (XSTAR  -  X(I) ) 

DIFF  ■  ABS(X (I)  -  XOLD) 

IF  (DIFF  .GE.  ERRMAX)  THEN 
ERRMAX  »  DIFF 
II  »  I 
END  IF 

SUM  *  SUM  ♦  DIFF 
110  CONTINUE 

C*********  CHECK  FOR  ITERATION  CONVERGENCE  ********************* 
ERRAVE  =  SUM/IMAX 

IF  (ERRAVE  .GE.  GRID ACC)  GO  TO  999 
200  WRITE (6, SO) K, ERRMAX ,11, ERRAVE 

56  FORMAT (/Tl 5, ’SOR  GRID  CALCULATION  ’, /,T5, ’USING  A  2ND  ORDER’, 

#’  UPWIND  DIFFERENCE  FOR  THE  CONVECTIVE  TERM’) 

80  FORMAT (T5, ’ ITERATION  #’ , I5,5X, ’MAX  ERROR  ■’,F10.3,’  AT’, 
#14,/,  SX, ’AVERAGE  ERROR  -’,F10.5) 

RETURN 

END 


SUBROUTINE  NEWP (THETA, C,NP,G) 

C*************************************************************** 

C  THIS  SUBROUTINE  CACULATES  THE  NEW  GRID  GENERATION  CONTROL  * 
C  FUNCTION  P.  IT  IS  BASED  ON  A  GLOBAL  TRUNCATION  ERROR  * 

C  ANALYSIS  FOR  THE  FLOW  SOLUTION  IN  THE  TRANSFORMED  PLANE  * 

C *************************************************************** 

COMMON  /A/IMAX, U(31 ) , ZI (51) , X (51 ) ,P(51) ,DX (51) , RE 
REAL  C<3),ZZ<11),UU(11),G(51) 

NN*(NP-i) /2 
WRITE (6, 13) 
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DQ  7  J-l.NP 
22  <J>  «ZI <J) 

UU(J)»U(J) 

7  CONTINUE 

CALL  ULSTSGXZZ.  UU, FLOAT <NP) , C) 

A2"C<3) 

Al-C<2) 

AO-C(l) 

CC  WRITE<6,27) <ZZ<J) , J-l.NP) 

CC  WRITE<6,28) <UU<J),J-1,NP> 

CC  WRITE (A, 29)  A2,A1,A0 

DO  S  I-2,NN+1 

G(I)  ■  A2*I**2  ♦  A1*I  +  AO 
DU»0.5t<U<I+l)-U<I-l)) 

DUU-U  <I+1)-2.0*U(I)+U(I-1) 

DUSIGN-SIGN ( 1 . 0, DU) 

IF(ABS(DU)  .LE.  .0000001)  DU-. 0000001 *DUSIGN 
DEN-2.0*A2II  +  A1 
DENSN-SIGN (1.0. DEN) 

IF  (ABS(DEN)  .LE.  0.0001  )  DEN-. 0001IDENSN 
POLD-P  < I ) 

P < I ) -POLD  +  THETAIDUU/DU  ~  THETAH2. 0HA2/DEN 
PMAX  -  U(I)*RE*DX(I) 

IF  <P<I)  .GT.  PMAX)  P ( I ) -PMAX 
IF  <P(I)  .GT.  2.0)  P(I)-2.0 
IF  <P(I)  .LT.  -2.0)  P(I)— 2.0 
WRITE  <6, 16) I , DU, DUU, DUU/DU, POLD,  P  < I > 

5  CONTINUE 

DO  10  I-2+NN, IMAX-NN-1 
DO  20  J-lpNP 

ZZ(J)  -  ZI (I+J-NN-1) 

UU<J)  -  U( I+J-NN-1) 

20  CONTINUE 

CALL  ULSTSQ <ZZ,UU, FLOAT <NP),C) 

A2-C(3) 

Al-C<2) 

A0-C<1) 

CC  WRITE (6, 27) (ZZ<J> , J-1,NP) 

CC  WRITE(6,28> (UU<J> , J-i,NP) 

CC  WRITE <6, 29)  A2,A1,A0 

G(I)  -  A2*H*2  +  Alii  +  AO 
DU-0.St(U<I+l)-U<I-l) ) 
DUU-Ua+l)-2.0*U<I)+U<I-l> 

DUSIGN-SIGN (1.0, DU) 

IF(ABS(DU)  .LE.  .0000001)  DU-. 0000001* DUS I GN 
DEN-2.0IA2II  +  A1 
DENSN-SIGN <1.0, DEN) 

IF  (ABS(DEN)  .LE.  0.0001  )  DEN-.0001*DENSN 
POLD-P ( I ) 

P ( I ) -POLD  +  THETAtDUU/DU  -  THETA*2. 0*A2/DEN 
PMAX  -  U<I) *RE*DX <I) 
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IF  (P(I>  .GT.  PMAX)  P ( I ) =PMAX 
IF  <P<I)  .GT.  2.0)  P(I)-2.0 
IF  <P ( I )  .LT.  -2.0)  P(I)»- 2.0 
WRITE (6, 16) I,DU,DUU,DUU/DU,PQLD,P(I) 

10  CONTINUE 

DO  30  J-NP,1,-1 

ZZ < J)  ■  ZKIMAX+J-NP) 

UU(J)  -  U<IMAX+J-NP) 

30  CONTINUE 

CALL  ULSTSQ ( ZZ , UU, FLOAT <NP> ,  C> 

A2«C<3) 

Al-C (2) 

AO-CU) 

CC  WRITE<6,27) (ZZ <J) , J»1 , NP> 

CC  WRITE (6,28) (UU(J) , J«1,NP) 

CC  WRITE (6, 29)  A2,A1,A0 

DO  25  I*IMAX-NN, IMAX-1 

G(I)  »  A2*I**2  +  A1HI  +  AO 
DU»0.5t(U<I+l)-U<I-l)) 

DUU-U(I+1)-2.0*U<I>+U<I-1) 

DUSIGN*SIGN ( 1 . 0, DU) 

IF (ABS(DU)  .LE.  .0000001)  DU«.0000001*DUSIGN 

DEN-2. 0*A2*I  +  A1 

DENSN-SIGN(1.0,DEN) 

IF  (ABS(OEN)  .LE.  0.0001  )  DEN-.OOOltDENSN 
POLD-P(I) 

P<I)«POLD  ♦  THETAIDUU/DU  -  THETAI2. 0* A2/DEN 
PMAX  «  U(I)*RE*DX(I) 

IF  <P(I)  .GT.  PMAX)  P(I)«PMAX 
IF  (P(I)  .GT.  2.0)  P(I)-2.0 
IF  <P(I)  .LT.  -2.0)  P(I)«~ 2.0 
WRITE (6, 16) I , DU, OUU, DUU/DU, POLD, P (I > 

25  CONTINUE 
P<1)»0.0 
G(1KI(1) 

G<IMAX)»U(IMAX) 

P(IMAX)»0.0 

15  FORMAT (/,T25, ’»***»***»  NEW  P  CALCULATION  ******’,// 
#T4,  ’ ZI’ , T14, ’ DU’ , T26, ’ DUU’ , T36, ’ DUU/DU’ , 

#T48, ’ POLD’ , T60, ’ PNEW’ , / ) 

16  F0RMAT(I5,8(3X,F9.3)> 

27  FORMAT (T2,’ZZ(J) i  ’ , 6F10.5) 

28  FORMAT <  T2,’UU(J)t  ’,6F10.5> 

29  FORMAT (T2/LST  SO  FIT  UU(ZZ)  » 

MF8.5, ’  *  ZII42  +  ’ ,F8.5, ’  *  ZI  +  ’,F8.5,) 

RETURN 
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